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In the vast majority of many-body problems, it is the kinetic energy part of the Hamiltonian that is best known 
microscopically, and it is the detailed form of the interactions between the particles, the potential energy term, 
that is harder to determine from first principles. An example is the case of high temperature superconductors: 
while a tight-binding model captures the kinetic term, it is not clear that there is superconductivity with only 
an onsite repulsion and, thus, that the problem is accurately described by the Hubbard model alone. Here we 
pose the question of whether, once the kinetic energy is fixed, a candidate ground state is groundstatable or 
not. The easiness to answer this question is strongly related to the presence or the absence of a sign problem 
in the system. When groundstatability is satisfied, it is simple to obtain the potential energy that will lead 
to such a ground state. As a concrete case study, we apply these ideas to different fermionic wavefunctions 
with superconductive or spin-density wave correlations and we also study the influence of Jastrow factors. The 
kinetic energy considered is a simple next nearest neighbor hopping term. 

PACS numbers: 71.10.Fd, 74.20.Fg, 02.70.Tt 



I. INTRODUCTION 

The problem of finding the ground state of a many-particle 
Hamiltonian is, in general, a daunting task. The problem is 
most severe in the case of fermionic particles, where the in- 
famous fermion sign problem plagues solutions via numerical 
methods. In contrast, in bosonic systems, Monte Carlo meth- 
ods are rather efficient in simulating systems of reasonably 
large sizes. Certain methods, such as Density Matrix Renor- 
malization Group, avoid the sign problem, but are mainly re- 
stricted to Id or quasi-ld systems. 

In this paper we step back from the problem of determin- 
ing the ground state of a given many-body Hamiltonian, and 
instead pose the following question: fixing the kinetic en- 
ergy part of the Hamiltonian, can a given wavefunction be the 
ground state for some choice of potential energy? Put simply, 
we ask if the wavefunction is groundstatable. The question is 
trivial to answer for bosonic system, as we discuss below, but 
in the case of fermions it is much more difficult and subtle. 

To illustrate this idea, let us start with a very simple exam- 
ple: the case of a single spin 1/2 degree of freedom. Consider 
a Hamiltonian of the form: 

H = -(j, + T4, (1) 

where is the usual spin-flip operator and Va is a diagonal 
operator in the {| t): I i)} basis. Here, we will not specify 
Va and try to diagonalize H; instead , we will consider the 
wavefunction 

= / [(1 - a) IT) + (! + «) li)] (2) 

and ask what the condition is on a so that \^a) is the ground 
state of H, for some a proper choice of Va ■ To answer that, 
the first step is to make an eigenstate of H by imposing 
H\'^a) ~ (so I^E'q) is an eigenstate of H with energy zero). 



The expression of Va follows immediately and we can rewrite 
iJ in a matrix form as: 

H={'-,^ 1^ . (3) 

The two eigenvalues of this problem are Ai = and A2 — 
+ 1^- Now, it is easy to see that |^'q,) will be the ground 
state of H if and only if a < 1 {i.e. if the wavefunction ele- 
ments are all positive). We will say that \^a) is groundstat- 
able for a < 1. On the contrary, when a > 1, I^'q) is an ex- 
cited state of the problem and no longer groundstatable. At the 
boundary between the two cases, one component of the wave- 
function vanishes at a = 1. Then, the potential energy blows 
up and the eigenvalue A2 goes from +00 to —oo. Of course, 
the property of groundstatability for a given wavefunction is 
directly related to the kinetic energy operator we have con- 
sidered. Had we chosen a different operator, we would have 
reached a different conclusion on a. The point is that once 
this operator is fixed, the problem is uniquely defined. 

This approach to the single spin Hamiltonian can be ex- 
tended to a many-body problem, where the kinetic energy is 
often chosen to be a local hopping operator between nearby 
sites. Then, from the set of all possible many-body wave- 
functions, some are groundstatable and others are not. It is 
of crucial importance to establish in which category a given 
wavefunction belongs to, since it determines if this state is al- 
lowed in nature. In general, for a given kinetic energy term, 
the Hilbert space is broken down in regions in which the wave- 
function satisfies groundstatability; as we will see below, the 
level of complexity of the partitions of the space into such re- 
gions is closely related to the presence or the absence of a sign 
problem in the Hamiltonian. 

The paper is organized as follows. In section we de- 
fine the problem of groundstatability on a finite dimensional 
Hilbert space and we show how it can be solved on a par- 
ticular case where the Hamiltonian admits a product form. 
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Then, we present in section |III] the main part of our work, 
namely how one can build a Hamiltonian for which a given 
many-body fermionic wavefunction is the groundstate. The 
numerical procedure is also detailed in that section. Results 
are shown in section HV] We discuss the case of the wave- 
function for non-interacting fermions which allows us to illus- 
trate the loss of groundstatability in these systems. We then 
present the potentials obtained from mean-field solutions of 
the Hubbard model, BCS superconductors and spin-density 
waves (SDW). These results are in accordance with mean- 
field analysis. By considering additional Jastrow factors, we 
also examine partially-projected BCS wavefunctions relevant 
for the study of high-Tc superconductors. Finally, a more open 
problem, with a class of wavefunctions containing both super- 
conductivity and antiferromagnetism, is investigated. 



this case is that ^'c > 0, VC. This is the case of matrix 
Hamiltonians for bosonic systems, and the strictly negative 
or zero off-diagonal elements is related to the absence of a 
sign-problem in the Hamiltonian. The problem of the single 
spin 1/2 mentioned in the introduction falls in this category. 
Now, one does not have the luxury of the stochastic mapping 
to a problem with positive probabilities in general. If some 
off-diagonal elements of the Hamiltonian are non positive, we 
lack any general theorem to conclude on the groundstatabil- 
ity of the wavefunction. Sometimes, it is possible to find a 
gauge transformation to come back to the simpler case with all 
strictly negative or zero off-diagonal elements; this happens 
for some spin models, like the anti-ferromagnetic Heisenberg 
Hamiltonian, where the sign structure of the ground state is 
given by the Marshall rule—. However, in fermionic problems 
and some frustrated magnets, the problem remains unsolved. 



II. THE PROBLEM 



A. General considerations 



Let us consider a finite dimensional matrix example. Take a 
Hamiltonian matrix i?c.cs where C, C index the states in the 
(iff -dimensional Hilbert space, for example the spatial config- 
urations of fermions on a finite lattice. Suppose that the off- 
diagonal elements Hc^c are known, and one wants to de- 
termine if the vector (state) with components '^c can be the 
ground state if the matrix elements in the diagonal are prop- 
erly picked. There are two steps in the problem; the first is 
trivial, to make | ^P) an eigenstate, and the second is the prob- 
lem we pose, whether it can be the ground state. 

We start by determining the diagonal elements from the 
condition that j^P) is an eigenstate. For simplicity, we shift 
again the eigenvalue to zero, and solve for the dn vari- 
ables Hcc in the diagonal: 



Hcc ~ ^ Hcc — - ^ Hcc — , 



c c^c 
so the Hamiltonian can be written as 



(4) 



H 



1 



cc 



c^c 



^\C){C\ + ^\C'){C'\ (5) 
^c ^c 

-\C){C'\-\C'){C\ , 



which is a sum of projector operators acting on a 2- 
dimensional subspace of states C, C" (one can check that the 
operator within brackets squares to a multiple of itself). 

The problem of groundstatability is that it is not guaranteed, 
with the Hamiltonian Hcc now determined, that |\E') is the 
ground state, and not an excited state. What are the conditions 
on the vector components 'i'c for it to be the ground state? 
If the off-diagonal matrix elements are all non-positive, then 
one can make use of the Perron-Frobenius theorem and the 
well-known connection to stochastic dynamics', or alterna- 
tively cast the Hamiltonian as a sum of positive semi-definite 
projectors^. Basically, the condition of groundstatability in 



B. Hamiltonians with a separable form 

Before going to our main case of interest which is the 
fermionic Hamiltonian with neighboring site hopping for ki- 
netic energy term, we would like to present another class of 
models where the question of groundstatability can be com- 
pletely and analytically answered even in the absence of a 
Marshall-like rule. The problem of groundstatability is de- 
fined for a given choice of kinetic energy operator, and we 
will choose here the off-diagonal elements of the Hamilto- 
nian H^fj, ~ ±wcwc> to be separable into products of 
wc, wc' € K- Notice that these models are highly non-local 
problems. However, they are interesting because they display 
the fundamental difference between the presence/absence of a 
sign problem in the system, mainly, the parameter space can 
be separated into disconnected groundstatable regions in the 
(+) case, whereas in the (— ) case the groundstatable region 
is just made of a single block. 

Let us consider the two possibilities: H^J^^j, = +wcwc' 

or H^(j^(j, = —wcwc- Now, given a vector "i!c, we con- 
struct the diagonal elements according to Eq. (HI so that ^c is 
an eigenvector with eigenvalue zero. The Schrodinger equa- 
tion for any eigenstate j?/'^) = X^c^c 1*-^) ^i'-^ energy A 
reads: 



± 



^ Wcwc Ipc' 

c^c 




■^c 

? WcWc 




and it is then straightforward to show that all eigenvalues A 
are solutions of the equation: 



Wr 



X±{Ec wc'fc)wc/'^c 



±1. 



(7) 



A = is indeed, by construction, a solution. The state 1 4*) is 
the ground state if all other solutions of Eq. (|7]i are positive. 
It follows that solutions of /-(A) = —1 satisfy A > if 
all poles of the function /_ (A) are positive, and solutions of 
/+(A) = +1 satisfy A > if one and only one of the poles of 
the function /+(A) is negative (see Figure[T]). 
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FIG. 1: (color online) Solutions of the Schrodinger equation ([7} (cir- 
cles) for a three dimensional Hilbert space when A = is the ground 
state energy. Notice the positions of the poles with respect to zero in 
the two cases. 

Notice that in the (— ) case, to fix the poles of /_(A) to be 
positive, the signs of groundstatable vectors 'i'c are related to 
those of wc, and one can thus write '^c = sgn(wc) l^cl^ 
which is a simple example of a Marshall sign. In this case, the 
signs of the uip's can be gauged away from the Hamiltonian, 
bringing it to the form that satisfy the Perron-Frobenius the- 
orem: Hc^c Hc^c = -\wc\\wc\ = 
|*cl >0. 

The condition for groundstatability in the (+) case is richer. 
The condition that one and only one of the poles of the func- 
tion /+(A) is negative leads to dn distinct sectors in the 
Hilbert space, each sector corresponding to the choice of 
which of the dn poles is selected to be the negative one. 
More explicitly, the condition on the poles is equivalent, for 
'^c ^ 0, to: 

«'c*c/E«^c'«'c') > Ofor(7 (8) 

C" 

wc^'c/(Xlwc7'^'c') <0 forCT^C (9) 

C" 

Each of these inequalities split the Hilbert space into two 
pieces via a hyperplane, and the dn conditions lead to a sim- 
plex, and the choices of the C to dn such simplexes (see Fig.|2] 
for simple examples on 3 x 3 Hamiltonians). In most of the 
cases, the edges of the simplexes correspond to the vanishing 
of one of the 'fc in the wavefunction. Indeed, coming from a 
groundstatable region, if one component \l/c changes sign the 
inequalities ^ are violated. On the edge, '^c = 0. the associ- 
ated diagonal element Hcc is infinite [see Eq. ^] and one of 
the positive eigenvalues diverges towards +00 and reappears 
at —00. Another edge is given by the equation ^ wc'^c = 0. 
In this case, the whole hamiltonian is reduced to the projec- 
tor H = \w) {w\ (this corresponds to the diagonal line on Fig. 
121) and there are only two eigenvalues: and 1. From this 
point of view, boundaries between groundstatable and non- 
groundstatable either corresponds to ill-defined Hamiltonians 
or to highly degenerate problems, the first case being the most 




FIG. 2: (color online) Examples of the domains where a wavefunc- 
tion I'l') = (^'i, 3*2, ^'3) is groundstatable: we fix = 1, and 
show the regions in the (*I'i,*I'2) plane where the wavefunctions 
are groundstatable. (— ) Hj^J^, = —1 (Frobenius case) and (+) 
h1^2„, — +1 (non-Frobenius case), for C, C" — 1, 2, 3. 



common one. 

For a Hamiltonian with off-diagonal matrix elements that 
cannot be separated into a product as in the example above, the 
situation is more complicated. Nonetheless, one should ex- 
pect that the feature we encountered in the simple case study 
should remain: generically the set of groundstatable wave- 
functions within the Hilbert space is largely fragmented into 
regions. In the case of the example, there are order dn, the 
dimension of the Hilbert space, regions. This is to be con- 
trasted to the case where there is no fermion sign problem, 
where there is one single region. This fragmentation of the 
groundstatable set should be a generic feature of systems with 
fermionic sign problems. Moreover, we also expect the diag- 
onal part of the Hamiltonian to abruptly change form when it 
leaves a groundstatable region and to become singular at the 
boundary (the edges of the simplexes). These abrupt changes 
in the Hamiltonian can be used as telltales that the wavefunc- 
tion, as function of some parameter, exits a groundstatable 
region. 

It is important to point out that, within a given groundstat- 
able region, many phases of matter can exist. Order param- 
eters computed from a groundstatable wavefunction can be 
used to classify the phases. The groundstatable regions thus 
do not delimit single phases; what they do demarcate are the 
regions where a wavefunction can possibly correspond to a 
state of matter, for a fixed kinetic energy term in the Hamilto- 
nian. 



III. MANY-BODY FERMIONIC HAMILTONIANS 

Let us now turn into a more practical application, and show 
how one can implement the procedure of finding the potential 
energy term for which a given fermionic many -body state 1 4') 
is the ground state, given a kinetic energy. We will consider 
specifically the case where the kinetic energy comes from a 
tight-binding hopping term on a lattice, which is common in 
many strongly-correlated electronic problems. We consider 
the case of fermions on a square lattice, and detail below the 
numerical procedure used to evaluate the potential energy of 
the Hamiltonian. 
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Shifting the ground state energy to zero, the Hamiltonian 

H^^V^-tJ2 clcja + H.c. (10) 



<IJ> 



that we seek should satisfy 

Hvtl*) = 



e„ > 0, 



(11a) 
(lib) 



for all eigenstates |A„). The i,a label the site and the spin 
of the fermions, respectively. (Hereafter we set the energy 
scale t = 1.) The potential Vj, depends only on the fermionic 
occupation operators riia- and is uniquely determined by the 
condition Eq. dl la| ) provided j'l') is groundstatable. We will 
focus here on the general form of this potential, addressing the 
question of groundstatability for more specific cases. 

We treat this problem in the configuration basis {|C)}, 
where a basis element stands for a set of positions of 
the 2N fermions, say {Rj};=i,...,Ar for the up spins and 
{R^}m=i N for the down spins. The anticommutation re- 
lations between fermionic operators also require enumerating 
the fermions and keeping the same ordering for each config- 
uration. The action of the kinetic operator T on an ordered 
configuration C can be understood by introducing a config- 
uration C' such that Hcc = {C'\f\C) = (C"|C"), where 
C is another ordered configuration differing from C by the 
local hopping of a single electron. If configuration C is cor- 
rectly ordered, Hcc is equal to —1, otherwise it is equal to 
+1. Thus, to determine the sign of {C'\T\C) one has to con- 
sider the positions of all the fermions in C and C". This is the 
sign problem, which bedevils Monte Carlo simulations on the 
Hubbard model in dimensions higher than one. 

In this study, the fermionic many-body wavefunctions will 
take the form: 



l*) = ^E"^cdet(^c) \C), 



(12) 



where Jc is a Jastrow factor that depends on the fermion oc- 
cupation numbers in configuration C, and (pc ^s, a N x N 
matrix with elements that depend on the position of the par- 
ticles in configuration C, [4>c]im = <y5 (^j — R-m)' '^^^^ V 
a function characterizing the correlations between pairs of 
fermions^. This form includes the wavefunction for non- 
interacting fermions, BCS and spin density wave (SDW) 
wavefunctions and also partial projections of these statesi^. 
In the configuration basis, the potential reads: 



{C\V^\C) = Vc 



E^cc'^$^, (13) 



Jc dot (0c) 



supposing there are no vanishing determinants. Only neigh- 
boring configurations, defined so that C and C differ by the 
hopping of a single fermion, contribute to the sum. To eval- 
uate the sum, rather than determining Hcc' for each pair of 
configuration C and C", we can use the fact that dct(0c") = 



{C\C') X dct(0p,) to rewrite the potential Vc as a function 
of configurations C: 



vc=-j:' 



Jc' det((?!)g,) 
Jc det(0c) ' 



(14) 



where the primed sum is over configurations C that dif- 
fer from C by the hopping of a single electron. We then 
compute Vc by calculating all the ratios det(0p,)/ det{(j)(=,) 
which are easy to calculate since the matrices 0^, differ from 
cf)c by the modification of one row or one column, depend- 
ing if an up or down spin hopped, respectively - recall that 

[4'c]im = f (^j — R-m^ Notice that, by working directly 

with the positions Hj and R^, issues of orderings of config- 
urations disappear from the problem. 

To evaluate the operator V^i, we first compute Vc for a large 
number of configurations M. The configurations are chosen 
according to their weight \'i'c\'^ via a Metropolis algorithm, 
using the inverse update method for fermionic Monte Carlo—. 
Then, we search for the best two-body approximation to the 
interaction, neglecting three-body and higher order terms: 

= H + Oifiihjhk) + 
H = f + E + uY^ n,^-n,^ + ^Y. ^-o ^^^3 (15) 

where fii = fii-^ + hii. (Notice that for fixed particle num- 
ber the onsite potential can be written with Va = U, with a 
constant shift absorbed into E.) The coefficients Vij are eval- 
uated through a linear least square method (LLS) and the best 
approximated solution is the one minimizing the sum of the 
squared residuals S = ((i?* — HY). This quantity can be re- 
lated to the overlap 5 between the wavefunction jvP), ground 
state of H, and the ground state |^) of H, as follows. Us- 
ing perturbation theory on 5H — H — H^, we can write the 
(unormalized) ground state wavefunction of H as 



i*> + E I") 



En 



(16) 



where \n) and En are the eigenstates and eigenvalues of H. 
The norm of this state, {tp\ip), can be related to the squared 
residuals S = - /f)^!^'): 



(V'iv) = i + E — 



E^ 



< 1 



Et 



1 + S/E'l 



(17) 
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Overlapping IvE*) with the normalized state j^*) ~ lip) , I , 

V {■4'H') 

yields (up to second order in perturbation theory in SH) 



(18) 



Therefore the squared residuals S* = {'i>\{H<}, — H)'^\'^) which 
we obtain by approximating the potential energy V^, by a two- 
body interaction are a measure of the overlap between the 
ground states of H^j, and its two-body approximation H. The 
smaller S, the closer the overlap is to 1. One can bound the 
overlap between the two wavefunctions by noticing that, even 
if the system is gapless, the excitation energy Ei should be 
controlled by the finite size L of the system, and thus if S is 
found to be much smaller than Ef, the overlap will remain 
close 1. In estimating the overlap hereafter, we use the worst 
case scenario that the system is gapless, with a wavevector 
2tt/L for the lowest energy excitation. (In the computations 
of S below, we assume linearly dispersing modes with a ve- 
locity of order unity.) 

We computed the potential and the associated overlap for 
several wavefunctions. In each case, we considered a tilted 
lattice of size + 1 with odd L and periodic boundary con- 
ditions to avoid singularities of d-wave wavefunctions in re- 
ciprocal space^. 10000 sweeps are usually considered for 
equilibration. Then, up to A/^ = 80000 configurations are 
taken for the LLS method. Computations have been made 
for various system sizes (L ~ 13, 15, 17, 19) at half filling 
(N = 170, 226, 290, 362). Because one can always change 
the constant Ehy a shift in all the coefficients Vij , we add the 
additional constraint J2\i-j\>R ^ij = such that the last L/2 
coefficients average to zero. We checked numerically that the 
results do not depend on that specific choice of R. 

We would like to emphasize here the difference between 
our method and the traditional Variational Monte Carlo 
(VMC) method. Given a wavefunction j^I'), the VMC method 
provides an upper bound for the ground state energy of a 
Hamiltonian H, i.e., in VMC the Hamiltonian is given, and 
a wavefunction is the target. By varying the parameters con- 
tained in l^*), one can find the best choice which minimizes 
the energy and then compute other operator averages such as 
order parameters or correlations functions. However, it is also 
possible that the real ground state of the system is so differ- 
ent from 1 5') that it cannot be reached by a variation in the 
parameters. This systematic uncertainty is not encountered in 
our approach as the Hamiltonian H^i derives uniquely from 
l^"), i.e., we inverted the target to be the Hamiltonian and not 
the wavefunction. If j^*) is groundstatable, it is the ground 
state of Hxji by construction. 



IV. RESULTS 

A. Fermions in ID 

The potential obtained from the LLS expansion should give 
an adequate description of the Hamiltonian when interactions 



are predominantly two-body. Accordingly, when three-body 
and higher interactions are totally absent, it should reproduce 
the exact form of the Hamiltonian. Hence, to check the consis- 
tency of our method, we would like to begin with a fermionic 
system with only two-body interactions whose ground state 
is known exactly, and try to recover the Hamiltonian starting 
from the wavefunction. Unfortunately, we lack any exact re- 
sults in two dimensions. So, we will preliminarily step back 
to the one dimensional case where exact results are known, 
and investigate systems of N interacting fermions on a ring. 
An appropriate case is the Hamiltonian with an inverse-square 
potential: 



H 



2A(A- l)7r^ 

Z2 



E- 

i<j sir 



1 



(19) 



Depending on the value of A, the potential can be either at- 
tractive (A < 1) or repulsive (A > 1), A = 1 corresponding 
to the non-interacting case. The ground state wavefunction of 
this Hamiltonian has been found by Sutherland^ to be of the 
product form 



^{xi,...xn) 



nsin> 

j>k 



TT{Xj - Xk) 



(20) 



Let us now apply our procedure to the wavefunction ( l20b . 
The configuration basis is the set of positions of N spinless 
fermions on a ring of size L. For M configurations |C), we 
calculate Vc using Eq. (|4]i and then perform the LLS method. 
We do not have to worry about groundstatability here, since 
the wavefunction is always the groundstate for any A. The 
potential energy we obtain is presented in Fig. |3]for differ- 
ent values of A, X = 302 and N = 80000. We plot the set 
of linear coefficients {[/, Vij} as function of the distance be- 
tween sites |i — j|. As shown, we recover a potential falling 
algebraically with the distance. The potential is attractive for 
A < 1, repulsive for A > 1, and vanishes in between (for 
A = 1). Indeed, the precise dependence of the potential en- 
ergy on the distance is found in quantitative agreement with 
that in the Hamiltonian [19] with our method returning an ex- 
ponent7 = 2.01±0.1 for the power law decay shown in Fig. [3] 
(bottom panel). For each value of A studied, the overlap 5 is 
found to be larger than 99.99%. 



B. Fermions in 2D 

In two dimensions, things are less simple for two reasons: 
first, except for very particular cases, we do not know if a 
wavefunction will be groundstatable or if it will be an excited 
state. Secondly, there are in general three-body and higher or- 
der terms in the potential that makes the expansion ( fTST i not 
exact and that will result in a reduction of the overlap 5. A 
good way to proceed then is to start from wavefunctions we 
know a priori are groundstatable {i.e. by other means) and 
then adiabatically deform them. By continuity, the resulting 
wavefunctions should also be groundstatable unless a bound- 
ary is met. Moreover, as long as the deviation of the overlap 
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FIG. 3: (color online)Top: Two-body potentials Vij as function of 
the distance \i — j\ between sites for the Calogero-Sutherland wave- 
function. Bottom: log V as function of log sin(7r(a;i — Xj)/L) for 
different values of A. 

away from 5 = 1 is small, the expansion ( fTSl l should be rele- 
vant. How can we detect a boundary then? Of course, we do 
not have any analytical criteria such as the inequalities (|9]l here 
but we can make some basic assumptions guided by what we 
learned from the separable case. At a boundary, the diagonal 
part of the Hamiltonian is ill-defined so we expect some kind 
of singularity in the Hamiltonian (notice that even the simple 
spin 1/2 example in the introduction displayed such singular 
behavior as one crossed the boundary of groundstatability). 
The singular behavior signaling the boundary of a groundstat- 
able region can appear in the set {U, Vij} or in the overlap 
S. In a finite system, that means a strong dependance of the 
Hamiltonian with the system size. Note the difference with 
a phase transition, where it is the wavefunction which dis- 
plays singular behavior at a critical point as the Hamiltonian 
is smoothly varied; here the problem is inverted, as it is the 
Hamiltonian that is singular at the boundary of the ground- 
statable region as the wavefunction is smoothly varied. 



1. The deformed non-interacting wavefunction 

Our starting point will be the Guztwiller wavefunction for 
non-interacting electrons. It is defined by Jc" = 1 and 

</'k(a < 0) = 1 

¥'k(a > 0) = (21) 
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FIG. 4: (color online) Representation of the function ipx in k-space. 
The dot line corresponds to the Fermi sea of the non-interacting sys- 
tem. 

where Lp\^ is the Fourier transform of (p{r) and = 
—2t (cos kx + cos ky) — fiQ , fiQ being the chemical potential^. 
It is the ground state of the tight-binding model where the op- 
erator V reduces to a constant. In fact, in one dimension, one 
can check, using Eq. ( fT3] l and the Vandermonde determinant 
formula that Vc is a constant independent of C. It can also be 
checked numerically in two dimensions. This type of wave- 
function is very useful because by changing the shape of the 
Fermi sea, one can also generate a set of excited states of the 
tight-binding model. These are, by definitition, non ground- 
statable. Having at our disposal a groundstatable wavefunc- 
tion and a set of non-groundstatable wavefunctions, we can 
ask the question how do we go from one to another This can 
be studied by considering the deformed wavefunction \'^x) 
defined by; 

= 1 0<\k\< ki 
ipl- = l-x ki < |k| < kp 
ifil ^ X kF< |k| < k2 (22) 
<^k = ^2 < |k| < ^, 

kp being the Fermi momentum. A particular choice of fci and 
/c2 is represented in k-space on Fig. ID By varying x from to 
1, we start in the ground state of the tight binding model and 
end in an excited state. During the process, we necessarly lose 
groundstatability. 

We studied the potential obtained from as function of x 
with the LLS method (see Fig. |5]l. For x ^ the two-body 
approximation shows a fast-decaying potential. We focused 
on the evolution of the Hubbard term U of this potential. As 
X increases from 0, U becomes first more and more negative, 
then abruptly changes sign around a critical value xc ~ 0.55, 
becomes largely positive, and finally steps back to zero. In- 
creasing the system size, the turnaround of U around xc be- 
comes more and more brutal. This behavior is also noticed 
in the other coefficients Vij. Another interesting feature is 
observable through the evolution of the overlap 5 as function 
of X (see Fig. |6] left panel). The overlap exhibits a grow- 
ing drop around xc, indicating the presence of large 3-body 
and higher order terms in the expansion of V . In the thermo- 
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FIG. 5: (color online) Left: Evolution of the potential for the deformed wavefunction ([23] with L — 15 and |i — j| < 5 at different values of x. 
At large distances, the potential identically vanishes. Full lines denotes potentials for which the wavefunction is the groundstate. Dashed lines 
indicate non groundstatable wavef unctions. Right: evolution of the Hubbard term U as function of x for different system sizes. 



dynamic limit, the brutal change in the form of the potential 
should eventually lead to a singularity in the Hamiltonian as 
function of x. We interpret the significant change in the nature 
of the potential and the fast increase of the correction to the 
overlap with system size near xc as signatures of the bound- 
ary of the groundstatable region. Precisely, for x < xq the 
state 1 4* a;) is indeed the groundstate of the (attractive) Hamil- 
tonian we are constructing, and for a; > xc it is just an excited 
state of the (repulsive) Hamiltonian. 

One can extract additional information on what is happen- 
ing near xc by probing the fidelity of the wavefunctio n ' ' . 
The fidelity, in this context, is a measure of the overlap be- 
tween two adjacent states in parameter space; 



(23) 



which has been proposed as a useful quantity to expose phase 
transitions". Indeed, at a critical point, F displays a drop that 
increases with system size, because the two states l^i^x) and 
|^'x+(5.t) describe two different phases of matter In the par- 
ticular case of a level crossing (first-order quantum phase tran- 
sition), the critical point also corresponds to a loss of ground- 
statability. 

We computed the evolution of the fidelity for \'^x) with 
6x = 0.005 (Fig. |6] top right). The fidelity does not dis- 
play any drop around xc- This fact suggests that the point 
Xc cannot be interpreted as a critical point (includying a first 
order transition). Instead, the situation appears to be that it 
is the Hamiltonian itself that becomes singular at xc (indeed 
much similarly to the simple case of the single spin 5 = 1/2 
discussed in the introduction). 



We observe two drops of at .t = and a; = 1. To un- 
derstand this, we measured the superconducting BCS order 
parameter: 



\m 




(24) 



ForO < a:: < 1, the system develops superconductivity (Fig. |6] 
bottom right). Like the fidelity, the SC order parameter is un- 
able to detect the loss of groundstatability at xc', the wave- 
function is continuous (again, it is the derived Hamiltonian 
that is not) and thus the order parameter derived from this 
wavefunction is non-singular at xc- But with our analysis 
of 6, we now know that, for x > xc, \^x) is not the ground 
state of the Hamiltonian that we constructed and so we can- 
not conclude on the presence of superconducting order in the 
ground state. Finally, note that one could have chosen a dif- 
ferent parametrization for the function and a different final 
excited state. For example, we checked several choices of ki 
and k2 ■ The evolution of U and S have been found to be simi- 
lar, just with different values for xc- 

The analysis of the Guztwiller wavefunction gives us the 
basic steps to follow in order to determine if a wavefunction 
is groundstatable: start from a wavefunction that is known 
to be the ground state of a Hamiltonian with a given kinetic 
energy, then change continuously a parameter and observe 
whether there is some rapid evolution of the potential and of 
the overlap. If no such feature appears, then the potential ob- 
tained from the LLS is indeed the potential for which | ^P) is 
the ground state of the Hamiltonian . 



J 



BCS wavefunctions 



where 



Introducing BCS pairing correlations between fermions, 
one can consider a superconducting wavefunction with 



(25) 



A (cos kx — cos ky), 
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FIG. 6: (color online) Left: Evolution of 5 as function of x for the 
deformed Gutzwiller wavefunction and different system sizes. Top- 
Right: evolution of the fidelity F as function of x for L = 15. 
Bottom-Right: evolution of the superconducting order parameter 
I < $ > I as function of x 



and (as in the Gutzwiller case) = —2t (cos + cos ky) — 
jjLQ. Here A and /^t are parameters (related, but not equal to 
the actual gap and chemical potential of the system). At half 
filling, we take ji ^ {) and vary the parameter A. We then 
compute {[/, Vij} and 5 for various system sizes. Let us first 
discuss groundstatability in this case. The BCS wavefunctions 
can be obtained adiabatically from the Guztwiller wavefunc- 
tion starting from A ^ 0. By adiabatically we mean that we 
did not find any singularity in either the Hamiltonian extracted 
from the wavefunction or the overlap 5 going from this limit 
to a finite A. This is understandable as we expect to open 
a gap by increasing A. Starting from a known groundstate 
wavefunction, we should remain groundstatable as long as we 
do not close the gap. Moreover, in the case of the BCS wave- 
functions, groundstatability is further supported by the fact 
that these wavefunctions are the ground states of a mean-field 
effective Hamiltonian. 

The potentials obtained for the s-wave state with L = 15 
are presented in Fig.|7]top. It shows a short distance two-body 
negative interaction whose strength is rapidly increasing with 
A. The potential vanishes when the fermions are separated 
by at least three lattice sites. So we see that the s wavefunc- 
tion seems a rather good approximation for the attractive Hub- 
bard model. The evolution of the overlap as function of A is 
shown in the Inset. We find 5 > 0.994 for A < 2.0, so the 
two-body approximation seems reasonable for these values. 
Surprisingly, the overlap seems to converge to a finite value 
as function of A(!). The two-body potentials for the d-wave 
case are presented in Fig. |7] bottom. The potentials show a 
complicate behavior as function of the distance with positive 
and negative coefficients in the limit of large A. The main 
common feature is the presence of a large negative nearest 
neighbor interaction term Vi, which is consistent with previ- 
ous mean-field analysis^^. Additional terms on a longer range 
are also non-zero due to the symmetry of the wavefunction. 
Comparing the magnitude of the potentials in the two cases, 
we find that the d-wave potential is always a lot weaker than 
the s-wave case. However, the overlap is smaller in the d- 
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FIG. 7: (color online) Two-body potentials Vij as function of the dis- 
tance jz — j\ between sites, evaluated at half filling. Top: the s-wave 
superconductor. Bottom: d-wave superconductor. Insets: The over- 
lap is always larger than 98.5% for the wavefunctions considered. 



wave than in the s wave case which means that the two-body 
approximation is less relevant for this symmetry. 



3. Partially projected BCS wavefunctions 

We also consider the partial Gutzwiller projections of BCS 
wavefunctions. These functions are defined (for the d-wave 
case) by (/3k = 'ft^ and a Jastrow factor 



(26) 



For < a < 1 this factor both penalizes double occupancy 
and is positive. Recalling the groundstatability conditions ^ 
obtained from the separable case, we expect that the wave- 
function remains groundstatable as long as the signs in the 
wavefunction are not changed. Notice also from (HI that a di- 
agonal element of the Hamiltonian becomes ill-defined (cross- 
ing from ±oo to Too) if the wavefunction ^'c for a config- 
uration C changes sign. Thus, the groundstatability of this 
projected BCS wavefunction is expected from by the fact that 
multiplication by a positive Jastrow factor does not change 
any signs of the BCS wavefunction. 

Figure [8] presents the evolution of the potential Vij for dif- 
ferent values of a and A = 0.5. The potential shows a large 
positive on-site potential growing with a and a smaller neg- 
ative short range interaction also growing with a. Thus, for 
values ofU^ 10 relevant for high-Tc superconductivity, we 
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FIG. 8: (color online) Potential Vij as function of distance for par- 
tially projected BCS wave function with Abcs ~ 0.5. Inset: evolu- 
tion of the overlap as function of a 



see that a projected BCS wavefunction is favored by a nearest 
neighbor attraction of order Vi ^ —3. A Hubbard term alone 
is not enough. Moreover, from the evolution of the overlap, 
we see that as a gets closer to 1, the two-body approximation 
is less and less justified. In fact, for a ~ 1 the wavefunction is 
the resonating valence bond state proposed by Anderson A 
better model to describe such a state would he at — J model 
where interactions are mediated through spin exchange. This 
is not allowed in our study since the Heisenberg part of the 
t — J model also contains off-diagonal interactions. The point 
a = 1 is a boundary where the partially projected wavefunc- 
tion loses its groundstatability. 



4. SDW wavefunction 

A spin-density wave state can be generally found in pres- 
ence of repulsive interactions, as expected from a mean-field 
solution of the Hubbard model. The SDW wavefunction is 
defined by: 



-*(Rj, R,i) = J2 ^(-a) «k(Rj) a_k(R|) (27) 



where the sum is restricted to the non interacting Fermi sea 
and 



with Q = (tt, tt) and: 



i(k+Q)-Rf 



(28) 




Is the SDW groundstatable? Yes, since we can recover the 
Gutzwiller limit by letting Ag^,^ without encounter- 
ing any singularity. The SDW wavefunction should remain 
groundstatable for any finite Ag^,^. The potential is pre- 
sented in Fig. |9l It shows a repulsive potential with a very 




FIG. 9: (color online) Two-body potentials Vij as function of the 
distance \i — j\ between sites, evaluated at half filling for the spin 
density wave antiferromagnet. Inset: The overlap is always larger 
than 99% The second inset for the SDW wavefunction shows the 
best fit for i = 25 and Agj^^ = 0.2. 



slow decay. The magnitude of the potential grows with Agj^^ . 
We find that the potential is well fitted by an algebraic decay 
V{r) ^ with /3 decreasing when Ag^^^ increases. We 
find 0.3 < (3 < 0.6 for 0.1 < Ag^.^ < 2.0. However, it 
is very hard to calculate this exponent with accuracy because 
our system is not large enough for the potential to really de- 
crease to zero. We also cannot totally exclude the possibility 
of a very large but finite range for the potential. The best fit 
for a system of size L ^ 25 with N = 313 and Agj^^ = 0.2 
is V{r) = exp(— r/5.9)/r'' '' (see second inset of figure|9l). 
Similarly to what we found in the case of the BCS wavefunc- 
tions, the two-body form is a good approximation for the po- 
tentials obtained from SDW wavefunctions as long as Ag^„ 
is not too large (the overlap 6 > 99% for Agj^^ < 2). 

Now, the method is not restricted to these simple cases. In- 
deed, It can be applied to any parametrization of (p{r). For 
example, one could search for the Hamiltonian for which a 
long-range wavefunction (e.g., with Lp{r) ^ l/r") is the 
ground state. The problem here will be not so much find- 
ing the Hamiltonian but knowing if we are actually starting 
the procedure from the real groundstate or from an excited 
state. Again, we should rely on some adiabaticity argument to 
answer this question. 



5. Mixed BCS-SDW wavefunction 

Several possibilities exist to construct wavefunctions with 
both BCS and AF orde r'"^''^-'^ . These in general rely on mean- 
field solutions of Hamiltonians having BCS and SDW cou- 
plings. Here, we will consider a different wavefunction de- 
fined by: 



R/,Ri„). (29) 



and we will take A^^g Ag^^ ~ 0.5. Although it is not ob- 
vious at first, this form also admits a decomposition in terms 
of single particle wavefunctions, as shown in the appendix. 
We would like to study the groundstatability of this wavefunc- 
tion and the evolution of the potential \4 as function of x. For 
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FIG. 10: (color online) Evolution of the potential as function of x around the pure SDW state 



very large positive or negative x, the wavefunction reduces to 
the usual s-BCS wave function dZSl l. For x = 0, this is the 
pure SDW wavefunction. The potential we obtain for an ar- 
bitrary X is presented in Fig.[TO] Starting from large negative 
values of x, the potential is attractive on a short distance. It 
does not vary much as soon as x < —0.3. When x is ap- 
proaching zero, the potential becomes more and more repul- 
sive. The transition is not smooth (see Fig. [TO]left). At some 
point, the Hamiltonian displays an attractive long-range po- 
tential with a short range repulsion around x = —0.06. This 
is quite unexpected since the s wave BCS corresponds to a 
short range attractive potential and the SDW to a long-range 
repulsion. Then, for a short range of very small and negative 
values of x (—0.03 < x < 0), the potential becomes purely 
repulsive with a Hubbard term larger than in the pure SDW 
case, i.e. U{x ^ -0.02) > U(x 0). Finally, from x = 
to X large and positive, the potential turns from repulsive to 
attractive in a very smooth way (see Fig. [TOl right). 

The evolution of the overlap S is shown in Fig. [TT] left. It 
displays a large drop in the whole region —0.2 < x < 0. But 



it is large and constant for all positive values of x up to zero. 
Another interesting information is given by the fidelity, shown 
in Fig. [TT] right. The fidelity is very close to 1 for all positive 
X. In contrast, it displays a sharp drop around x = —0.06 
which grows with the system size. 

The fact that nothing happens in both the fidelity and the 
overlap 5 for x > leads us to the conclusion that the process 
of going from x = +oo to .t = preserves the groundstata- 
bility of the wavefunction. Therefore, we are able to find an 
adiabatic path between the BCS state and the SDW state. In 
contrast, for x = — oo to .t = 0, we face a phase transition 
near x = —0.06; this transition is probably first-order, given 
1) the sharpness in the drop of the fidelity, and 2) the fact that 
the states at x = ±00 are the same (they differ by an overall 
sign of the wavefunction) and thus cannot be separated by a 
second order transition. This particular evolution is peculiar 
to the mixtured considered in ( |29l ). For instance, the wave- 
function proposed by Giamarchi and Lhuillier in Ref.[T4ldoes 
not display this behavior. 



r 



To gain some insights from what happens close to the phase 
transition, we also measured the two order parameters; the AF 
order parameter m defined by 



-E 



(30) 



and the SC order parameter (l24l) . The evolution of the order 
parameters as function of x are presented in Figure [T2] Let 
us first discuss the evolution of the AF order parameter. Start- 
ing from X large and negative, the magnetization steadily in- 
creases from zero. It then displays a maximum at x = —0.06 
and finally decreases back to zero for x large and positive. So 
we find that the maximum of the magnetization does not cor- 
respond to the pure SDW state but rather to the SDW state 
with small additional BCS correlations. The study of the 
superconducting order is also very interesting: ($) is maxi- 



mum for large values of |x| and it vanishes at x = 0.0 as ex- 
pected. But It also displays an unexpected local maximum at 
X = —0.06. In the range —0.12 < x < 0, both superconduc- 
tivity and magnetism are not competing but rather supporting 
each other. Notice also that the state with x = —0.12 and 
the pure SDW state share the same characteristics: they have 
the same value of to and no BCS order at all. What does this 
corresponds to in terms of the Hamiltonian? Turning back to 
Fig-Sni we see that the maximum of both orders corresponds 
to a potential with a short range repulsion and a long-range 
attraction. However, It is hard to draw a clear conclusion on 
this potential since the variation of the overlap 8 is pretty large 
near x = —0.06. Nonetheless, it appears that we can trust the 
results on the potential for —0.04 < x < where the over- 
lap is still large. As we already discussed, this corresponds to 
a purely repulsive potential with a Hubbard term larger than 
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FIG. 12: (color online) Evolution of the staggered magnetization per 
site and the SC order parameter for L = 13 as function of x. Inset: 
Evolution of the Fidelity as function of x 



in the pure SDW case. A more detailed investigation of this 
phase could give interesting results on the possibility of hav- 
ing BCS order with purely repulsive Hamiltonians. 



V. CONCLUSION 

In conclusion, although the question here posed of ground- 
statability of a given wavefunction for a fixed kinetic energy 
term is much harder to address in fermionic than in bosonic 
systems, there are cases where it can be answered concretely, 
and we gave examples in this paper. Considering exact re- 
sults such as the Calogero-Sutherland wavefunction in ID, 
we were able with our method to recover the Hamiltonian 
starting from the wavefunction. In 2D, we illustrated the 
problem of groundstatability on the Guztwiller wavefunction, 
where by slightly deforming the Fermi sea of a non-interacting 
fermion system, one can go from the groundstate to an excited 
state. Then, we analyzed several mean-field wavefunctions 
with different types of superconductivity and antiferromag- 
netism. These wavefunctions appeared to be groundstatable 
as they can be obtained by a deformation of the Guztwiller 
state without losing groundstatability in the process. Starting 
from these wavefunctions, we obtained the potential for which 
these states are the exact ground states, and we showed that the 
two-body approximation to the potential appeared to be valid 
in a broad range of parameter space. We found both short- 
range and long-range interactions in the Hamiltonian. Finally, 



we discussed the case of two non-trivial cases: the partially 
projected Guztwiller wavefunction and a state with a mixture 
of superconductivity and antiferromagnetism. In the latter, we 
were able to find a Hamiltonian favoring both SC and AF at 
the same time. 

The approach we follow, of constructing the Hamiltonian 
starting from the wavefunction and kinetic term, should be 
useful in determining whether certain exotic states of matter 
- for example non-Fermi liquids in two or higher dimensions, 
and RVB states - are permitted in nature. Instead of guessing 
Hamiltonians that would realize these states, one can algo- 
rithmically determined the target Hamiltonian starting from a 
wavefunction and local kinetic terms. Whether such states of 
matter exist in nature translates into the question of whether 
these wavefunctions are groundstatable or not. 

We thank F. Alet and A. Sandvik for enlightening discus- 
sions, and GENCI for allocation of CPU time. Simulations 
used the ALPS librariesi^. 



APPENDIX A: BUILDING THE dot [(^] MANY-BODY 
WAVEFUNCTION FROM ONE PARTICLE 
WAVEFUNCTIONS 

It is possible to generalize the results in Refs4i^ and write 
many-body wavefunctions (built by creating particle pairs) 
in terms of a determinant of a matrix built from functions 
ip{R^ ,R'^) of two variables, the positions and of up 
and down particles. It is actually interesting to ask the reverse 
question, and find out the conditions on a function (p{R^ , R^) 
so that the determinant of a matrix constructed from this func- 
tion corresponds to a many-body wavefunction built by creat- 
ing particles in pairs. The reason for addressing this question 
is that one can then use such function (p{R^ , R^) to contruct 
interesting many-body states where different types of order 
co-exist. 

Consider the state 

where N is the number of pairs, and |0) is the empty (vacuum) 
state. The wavefunction is given by 

({Rj},{Ri}|M')=dct[0] (A2) 

where the N x N matrix [(j)\im = ^ {^\-, R-m) ' ^^'^ the func- 
tion (f is given in terms of the single particle wavefunctions of 
the states created by c\ ^ and Cj(._^^ ^ and labeled by A: 

^(RT,R^) = E 0a al(RT) (R^) . (A3) 

A 

Now, let us suppose that we want to start with a func- 
tion tfirO , R^), and determine how to decompose it in terms 
of single particle wavefunctions q:]^(RT) and aj.(.^j(R^) as 
above. We will do this construction as follows. 
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The positions and R^^ take values over the Ng lattice 
sites Ra, a = 1, . . . , Ng. We can thus consider the function ip 
as a Ng X A^^ matrix (pab = fCR^ = R-a, = Rfe)- First, 
note that the matrices ip'^ip and ipip^ are hermitian and thus 
diagonalizable: 

ip(p'' al = ex al (A4) 

That the indices A labeling the states and the eigenvalue are 
common can be seen as follows: 

((^^t) Q,T ^ Q,T ^^5) 

and therefore ((^^ q,T ~) jg eigenstate of <^^^<^ with eigenvalue 
e\. Indeed, we can actually pair up the eigenvalues of ip^ip 
and (fiip^'- (x {(p'^ a\) and a\ oc {(p a\j. More precisely, 
we can write {(p'' a\) = and a^) = (f'xctx^ where 

I'/'Ap = £a- (Notice that the phase of (px can be removed by 
choosing the overall phase of the eigenvectors.) 
We can thus construct an operator 

^ = ^0A , (A8) 

A 

from which we can write back the function (^(R^jR-^) by 
sandwiching (^between (R^| and |R^): 

<p(RT,R^) = 5] 0a ai(R^) ai*{R^) . (A9) 

A 

What we now need is a symmetry a| (R'^) = '^/(a)(-^'^) 
that enables us to identify Eq. ( IA3l l with ( IA9b . 

Notice that if ipiR^ , R^) is a real function, and if the spec- 
trum of p^p is non degenerate, then its eigenvectors q;| are 



necessarily real and one has (R^) = ce\{'R^)- So to get 
a non trivial function /, it is important that all eigenvalues of 
p}^Lp are degenerate (except maybe at some particular point 
where /(A) = A), and the corresponding eigenvectors reside 
in an eigenspace E\ of dimension 2. For real p it is possi- 
ble to write two orthonormal real eigenvectors in E\, which 
can then be used as real and imaginary parts of new orthogo- 
nal eigenvectors and . We then identify the conjugate 
state = o,\:(^\y where the pair A, / (A) labels the two states 
in Ex = Efi^x)^ completing the construction. 

What do we need to make the spectrum of p^p degener- 
ate? Suppose the function (^(R^,R^) has some symmetry. 
For instance, in the case of the BCS and SDW wavefunc- 
tions, one can check that (^(R1^,R-^) = <^(— R^,— R^-) is 
indeed a symmetry. This symmetry operation, in terms of 
the matrix p, is implemented through a Hermitian operator 
P such that PpP ~ p, or equivalently [P, p] = 0. It can be 
trivially checked that \P, p^ ~ as well, and consequently 
[P, (^t^] = 0. Because P has two different eigenvalues ±1, 
the eigenspaces E\ have dimension 2, which is exactly what 
we need to construct the corresponding / (A) to a given A with 

I I * 

Qft/xN = at . Therefore, under the conditions above, the 

iW X 

identification of Eq. ( |A3t with iA9i is complete. 

Finally, notice that if one assembles a function from a 
linear combination of two functions that satisfy the condi- 
tions above (for example, a symmetry such as pilO, R^) = 
(^(-RT, -RJ^) as in the BCS and SDW cases), then the re- 
sulting function also satisfy the conditions. In particular, the 
combination 

P^{Rj,Rl) = X pl^.iR] ~ RL) + ^sow(R/,Ri)- 

(AlO) 

does lead to a good fermionic wavefunction (built as in 
Eq. ( IA2b ); this type of wavefunction is the starting point to 
the studies that we carried in section llVB5l 



* C. Castelnovo, C. Chamon, C. Mudry and P. Pujol, 
Ann. Phys. (N.Y.) 318, 316 (2005). 

^ D. P. Arovas and S. M. Girvin, in Recent Progress in Many Body 
Theories, edited by T. L. Ainsworth et al. (Plenum, New York, 
1992), Vol. 3, pp. 315-344. 

^ W. Marshall, Proc. Roy. Soc. (London) A 232, 48 (1955). 

" C. Gros, Phys. Rev. B, 38, 931 (1988). 

^ J.P. Bouchaud, A. Georges and C. Lhuillier, J. Physique, 49, 553 
(1988). 

D. Ceperley, G.V. Chester and M.H. Kalos, Phys. Rev. B, 16, 3081 
(1977). 

' A. Paramekanti, M. Randeria and N. Trivedi, Phys. Rev. B 70, 
054504 (2004); C. Gros, R. Joynt and T. M. Rice, Z. Phys. B 68, 
425 (1987) and Phys. Rev. B 36, 381 (1987). 

* B. Sutherland, Phys. Rev. A, 4, 2019 (1971) 

' M. C. Gutzwiller, Phys. Rev. Lett 10, 159 (1963); Phys. Rev. 137, 



A1726 (1965). 

W. L. You, Y. W. Li, and S. J. Gu, Phys. Rev. E 76, 022101 (2007). 
" P Zanardi, M. Cozzini and P Giorda, J. Stat. Mech. (2007) 
L02002 

'2 E. Dagotto et al., Phys. Rev. B, 49, 3548 (1994). 
PW. Anderson, Science 235 1196 (1987). 
T. Giamarchi and C. Lhuillier. Phys. Rev. B, 43, 12943,(1991) 
L. Spanu, M. Lugas, F. Becca and S. Sorella, Phys. Rev. B, 77, 
024510 (2008) 

Cedric Weber, Andreas Laeuchli, Frederic Mila, Thierry Gia- 
marchi, Phys. Rev. B 73, 014519 (2006) 
" R Albuquerque et al., J. Magn. Magn. Mater. 310, 1187 (2007); 
M. Troyer, B. Ammon and E. Heeb, Lect. Notes Comput. Sci., 
1505, 191 (1998). See ihttp : //alps . comp-phys ■ oTjl 



